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Abstract 

Laser speckle has been proposed in a number of papers as a high-entropy source of unpre- 
dictable bits for use in security applications. Bit strings derived from speckle can be used for a 
^ , variety of security purposes such as identification, authentication, anti- counterfeiting, secure 

key storage, random number generation and tamper protection. The choice of laser speckle as 
a source of random keys is quite natural, given the chaotic properties of speckle. However, this 
same chaotic behaviour also causes reproducibility problems. Cryptographic protocols require 
either zero noise or very low noise in their inputs; hence the issue of error rates is critical to 
applications of laser speckle in cryptography. Most of the literature uses an error reduction 
method based on Gabor filtering. Though the method is successful, it has not been thoroughly 
analysed. 

' In this paper we present a statistical analysis of Gabor- filtered speckle patterns. We in- 

^ ■ troduce a model in which perturbations are described as random phase changes in the source 

! plane. Using this model we compute the second and fourth order statistics of Gabor coef- 

VO ] ficients. We determine the mutual information between perturbed and unperturbed Gabor 

coefficients and the bit error rate in the derived bit string. The mutual information provides 
^ . an absolute upper bound on the number of secure bits that can be reproducibly extracted from 

I noisy measurements. 

Keywords: Physical Unclonable Function, PUF, speckle, Gabor transform, entropy, key 
extraction, fuzzy extractor 
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1 Introduction 

1.1 Key generation from speckle patterns 

In [HI [12] Pappu et al. proposed to use speckle patterns obtained from coherent multiple 
scattering in a token to authenticate persons and devices. In a typical scenario, a person 
carries an authentication token consisting of a transparent material with scattering particles 
inside, e.g. glass with air bubbles. When he wants to get access to some service, he presents 
his token to a reader device. The device shines laser light onto the token under some pre- 
determined conditions (wave length, angle, focal distance, beam shape etc). This is called a 
'challenge'. The resulting speckle pattern (in transmission or reflection) under some prede- 
termined angle is recorded by the device. The recorded image is called the 'response'. The 
response is processed, e.g. by Gabor filtering, to yield a bit string that is reasonably insensi- 
tive to noise in the image. This bit string is compared to a previously enrolled bit string. If 
the strings are sufficiently similar, the token is authenticated. (A variant of this procedure, 
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involving a scan over a length of paper, was developed in [T] for the authentication of paper 
documents.) 

Often one cannot trust the reader device and/or the link between the reader and the 
verifier. In that case a Challenge-Response Pair (CRP) cannot be used safely more than once 
in the above scenario. In order to have a secure and practical token, it must be possible to 
obtain many different CRPs from one token. It must also be very hard to predict a CRP given 
previously observed CRPs. Further security requirements follow if one demands that it must 
be hard for an attacker to (i) extract all CRPs from a token in a short amount of time, and (ii) 
to extract enough information from the token either to physically clone it or to successfully 
compute its responses. Pappu introduced the name 'PUF' for a token (not necessarily optical) 
that satisfies all these security requirements. PUF stands for Unclonable Physical Function. 
Alternative names in the literature are Physical One- Way Function (POWF) and Physical 
Random Function. The word function stems from the fact that a response can be regarded 
as the evaluation of a complicated function of the argument; the function is parametrised by 
the physical structure of the token. It turns out that the physics of multiple scattering is 
compatible with all PUF requirements, especially if the token is created by a random mixing 
procedure of sufficiently small particles. 

Going one step further than the simple matching procedure of [HI [12] , it is possible to 
use a token's response as a secret key in a cryptographic protocol. This is nontrivial, since 
any amount of noise is fatal to ordinary cryptographic primitives. Secure forms of error 
correction, in which the redundancy data does not leak (much) information on the secret 
key, were developed in [9l H] . These techniques are called fuzzy extractors or helper data 
schemes. Their application to optical PUFs was studied in [17[ IT3l I14j . Key generation 
from CRP measurements is an enabler for a wide variety of security applications such as 
authentication, brand protection, tamper protection, anti-counterfeiting, secure key storage 
and special forms of authenticated computation [6j. For an overview of the subject of security 
with noisy data we refer to |15j . 

In all of these examples, it is important to have a good understanding of the number 
of random bits that can be extracted from the measurements. Overestimation can lead to 
serious cryptographic weaknesses. Underestimation leads to waste of resources. A general 
framework for the computation of measurement entropy was set up in [TB] and applied to 
transmissive optical PUFs. 

A different approach was taken in [8], where Gabor- filtered speckle patterns were com- 
pressed using the Context Tree Weighting (CTW) method. The size of the compressed data 
gives an upper bound on the entropy. 

A second important point is a good understanding of the noise that occurs in the response 
when the same challenge is applied multiple times. This noise determines how much of the 
total entropy of a response can be extracted in a reproducible way. Measurement noise 
is caused by many factors: temperature, moisture, stray light, mechanical misalignment, 
differences between reader devices, ageing etc. In [17J several methods were proposed to deal 
with noise in optical PUFs, e.g. alignment methods and efficient protocols. In [8] CTW 
compression was employed to estimate the mutual information between two (Gabor-filtered) 
noisy measurements of the same response. This information-theoretic quantity captures the 
shared entropy between two data sets and gives an upper limit on the length of the shared 
key that can be reproducibly extracted from these sets. 

That work has resulted in a lot of practical know-how, sufficient to set up a secure key 
extraction system. What is lacking, however, is a theoretical understanding of the effects 
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of measurement noise on the Gabor coefficients. In this paper we will address the issue of 
random perturbations and the statistical properties of Gabor-transformed speckle patterns. 

1.2 Binarized Gabor coefRcients 

Bit strings can be extracted effectively by using a Gabor transform. This method was pro- 
posed in [11] and further studied in [T71 [8]. Gabor Transforms are well suited since they 
are insensitive to small changes in an image and they reveal the locations as well as the 
orientations of structures at different spatial frequencies. They are used in a wide range of 
applications, such as iris recognition [3J, texture analysis and image enhancement, coding and 
compression. 

Here we briefly review the method used in the literature. A laser beam illuminates an 
object and, either by transmission or reflection, produces a speckle pattern. An image of the 
pattern is recorded in the 'detection plane'. A point in the detection plane is denoted as 
two-dimensional vector x. The light intensity in the detection plane is denoted as I{x). 

A two-dimensional Gabor basis function T{w, k, xq,x) is the product of a plane wave with 
wave vector k and a Gaussian with width w centered on xq. We write the Gabor basis 
functions Tjm and the Gabor coefficients G as follows: 



Only the imaginary (sine) part of the transform is considered. This is motivated by the 
property that the imaginary part is invariant under spatially constant perturbations of the 
intensity. 

Gabor coefficients G are evaluated for a subset of parameters w, k, xq, e.g. on a sub- 
lattice of positions xq, for two perpendicular choices of k (with equal modulus \k\), and for 
one fixed w. Since the basis functions form an over complete set, such a restricted choice of 
parameters can capture almost all information available in a speckle image. 

CoefRcients are discarded if they do not exceed a certain threshold T, i.e. one only keeps 
|G| > T. The chosen coefficients are called 'robust', because they are unlikely to be affected 
by noise. Finally, the robust coefficients are binarized; positive values are mapped to '1' and 
negative to '0'. 

This procedure is applied to the speckle pattern photographed during enrollment, and 
again to the image obtained in the authentication measurement. The Hamming distance 
(number of bit flips) between the enrolled bitstring and the second bitstring depends on the 
threshold T and the amount of measurement noise. Based on knowledge about the expected 
number of bit flips, one applies an error correction scheme that can cope with the noise. It 
is important to keep in mind that the robust bitstring can be deceptively long. The actual 
amount of information contained in it can be much less than the length, due to correlations 
between the Gabor coefficients [17j . 

1.3 Contributions and outline of this paper 

In Section [2] we first briefly describe the random phase model and the intensity statistics that 
are obtained from it. We motivate our use of this model. In Sections [3] and H] we study the 
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statistics of Gabor coefficients and the effects of random perturbations. This paper contains 
the following novel contributions: 

• In Section 13.21 we analyse the statistical properties of the Gabor-transform ([1]) of a 
speckle pattern. We present a procedure for computing arbitrary moments of the dis- 
tribution. Computation of the first four moments shows that there are small deviations 
from the normal form. 

• In Section 13.31 we compute the information content of a set of Gabor coefficients, for a 
given noise level of the detector. The entropy per typical speckle area turns out to be 
proportional to the square of the logarithm of the signal to noise ratio. 

• In Section 14.11 we introduce a method of perturbing a speckle pattern in the random 
phase model. Each sized source region has its phase shifted by a small random 
amount e, where e is drawn from a uniform distribution of width 2q. By tuning q £ 
[^Aip, vr] (where Aip denotes the phase uncertainty due to the number-phase uncertainty 
relation) , the magnitude of the perturbation is selected. We use this kind of perturbation 
to represent a misalignment, such as a shift or rotation of the token or the laser, or a 
change in the structure of a token. 

• In Section [4.21 we compute the mutual information between the original speckle 'source' 
and the perturbed one as a function of the noise strength q. The result is proportional 
to log(7r/g). 

• In Section 14.31 we show that the correlation between the intensity before and after a 
perturbation is given by ^^^r^- In Section [4.41 the same correlation is obtained between 
Gabor coefficients before and after a perturbation. 

• In Section [4.4.1l we compute the mutual information between a set of Gabor coefficients 
before and after a perturbation. 

• In Section 14.4.21 we compute the bit flip probability for a binarized Gabor coefficient 
due to a random perturbation. 

• In Section [5] we give experimental results, (i) It turns out that the empirical distribution 
function of the Gabor coefficients is consistent with theory. There is a noticeable devi- 
ation from the Gaussian form. The theoretical prediction of the variance matches very 
well with the data, (ii) We studied random perturbations by doing measurements on a 
sample whose surface structure slowly changes in time. Correlations were determined 
between the state before and after a perturbation. As expected from the theoretical 
results, there is a linear relation between the correlation function of the intensity and 
the correlation function of the Gabor coefficients. 



2 The random phase model 
2.1 Motivation and definitions 

Throughout this paper we use the random phase model as described by Goodman [7], in the 
Presnel approximation, in a free space geometry, for completely polarised light. This model 
has the advantage of being relatively simple while yielding intensity statistics that agree with 
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experimental observations. We depart from the traditional approach only in one respect. 
In [7] the components of the electric field amplitude {A^, Ay) in the detection plane are 
sometimes treated as the 'fundamental' degrees of freedom. For instance, second order and 
higher order intensity correlations are derived using the Gaussian distribution of A{x). The 
traditional approach has the drawback that it is very difficult to keep track of the number of 
degrees of freedom: It looks as if there is one degree of freedom per (continuum!) location x 
in the source plane. In fact the physical degrees of freedom lie in the source plane (defined 
as the exit plane of a transmissive PUF, or as the surface that reflects the laser light) and 
they are very easy to identify and to count. Another drawback is that there is no natural 
way to introduce misalignment perturbations in terms of the A{x) variables. On a more 
esthetic level, there is the drawback of having to rely on the Central Limit Theorem to get 
the Gaussian distribution of A, while in fact the number of random amplitudes added together 
is not infinite but merely very large. 

For these reasons we base our calculations on the random phases in the source plane as 
the fundamental degrees of freedom. All the well known speckle properties are of course 
reproduced in this approach. The model looks as follows. Diffused light leaves the PUF at 
the exit plane, through a disc-shaped region with radius R which we call the 'source'. We 
assume that the intensity is the same everywhere in the source. (We normalize the intensity 
to 1). Hence the source is modelled as a collection of random phases 99. The disc is divided 
into small regions of area A^. (The total number of regions is denoted as A'^reg = ttB? /}?). 
Together these generate the speckle pattern according to Huygens' principle. The complex 
amplitudes a in each region are the basic degrees of freedom. We write 

= exp (3) 

where the subscript a denotes a discrete two-dimensional coordinate in the source, with 
\a\ < R. The phases ipg at all the locations a are independent stochastic variables, with 
a uniform distribution in the interval (— vr, vr]. We introduce the notation (•)^ for taking the 
expectation value with respect to the random phases. We have 

{{a,r)^ = fornGN+ ; (a.^ai)^ = (4) 

Prom these basic rules it is straightforward to derive many-point correlations. 

Note that in adopting independent random phases we ignore the correlations that are 
known to exist between the phases as a consequence of either (a) multiple coherent scattering 
in a diffusive medium (see e.g. [5]), or (b) height and/or orientation correlations between 
microscopic pieces of a rough surface. These correlations lead to a reduction of the number of 
degrees of freedom. However, in this paper we are primarily interested in the influence that 
the parameter choices in ([2]) have on the information that can be extracted from a Gabor- 
filtered speckle pattern. In this context the correlations between phases in the source plane 
are only of minor importance. Hence we will ignore them and work with dH). 

The distance between the source and the detection plane is denoted as z. We assume 
z 3> A and use the Fresnel approximation. For the complex amplitude A = A^ + lAy in a 
point X in the detection plane we then have 

A{x) = — Qgexp —i—{x — g)^. (5) 
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2.2 Entropy of the source 



The entropy of the source is an upper bound on the entropy of the speckle pattern in a half- 
sphere. The entropy of the source is easily computed in the random phase model. The phase 
distribution is completely uniform. Therefore the entropy reduces to the logarithm of the 
number of possible states. It is well known that a coherent state with photon number Nq 
has an uncertainty AA^ = \/No in the photon number. Using the number-phase uncertainty 
relation ANAip = ^, we find a phase discretisation Aip = l/{2^/I%). Hence the entropy 
(expressed in bits) is given by 

077- 

H[a] = log2( — )^-^ = iVreg log2(47ryiVo). (6) 

Note that this result is equivalent to the estimate in [T7] , where the light exiting the PUF was 
described in terms of transversal momentum modes. In [T7] the correlations between modes 
were studied as well, and the resulting entropy reduction was estimated. In this paper we 
will not take such correlations into account. 

In order to get some feeling for the orders of magnitude we substitute numbers into ([6|) . A 
laser with A = 780nm produces an output power P =lmW, and a measurement with a CCD 
camera takes about At =lms. The total number of photons involved in one measurement 
is PAt/(hc/X), where h is Planck's constant and c is the velocity of light. Thus we arrive 
at A^o = XP At / (hcNj-Qg) . Assuming a source diameter of 1mm, we get iVreg = 1.3 • 10^ and 
No = 3- 10^ photons per region, yielding an entropy of approximately 14 bits per region of 
size A^. 



2.3 Statistics of the intensity 



All the well known statistical properties of the intensity can be derived from the random 
phase model. For completeness and for use in later sections, we briefly discuss how these 
properties are derived. The intensity at position x in the plane of detection is given by the 
squared modulus of the amplitude ([5]), 

I{x) = \A{x)\'^ = ~2 X] 



Xz 



W -a^ + 2x-{a-h) 



(7) 



a, 6 



The average intensity lav is obtained by taking the expectation value (•)^ and directly apply- 
ing Summation over a constant yields a factor A^reg- We obtain 



For second order statistics of the intensity one needs 4th order correlations of the random 
amplitudes a. In particular, from ^ it follows that 



a\,b\ 02,02 ai,02 12,01 



Using ([UJ, the well known results follow for the variance (cr/ = /av) 
correlation function C/, 



(9) 

and for the intensity 



Ci{x,x) :-- 



{mm). 



Ji{\x' -x\/M) 
\x' - x\/M 



-\ 2 



(10) 
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where M is a constant proportional to the average speckle size, 



M 



27rR 



(11) 



and Ji is a Bessel function. Higher order expectation values can also be computed. In 
particular, in order to derive the well known exponential probability density 



it has to be shown that 
function, which also follows from ([j 



p(/)=/-,iexp(-//4v), (12) 
/^n!. This is done using the following correlation 



Q, 



ai 



. a 



En*. 



ueSn k=l 



(13) 



Here Sn stands for the 'symmetric group' of all the n\ possible permutations of the numbers 
1 . . . n. Eq. (jlSp is also used in the derivation of the joint probability distribution p(/(x), I{x')). 
This distribution follows from the expectation value 



[/(f)r[/(x')r>^ = C+^'"n!m! 2F,i-n,-m;l;Cjix,x')) 



and is given by 

p(/(f),/(f')) 



11(1 -Ci) 



exp 



I{x) + I{x') 
' 4v(l - Ci) 



2^i{x)i{x') ycy 



(14) 



(15) 



Here C/ is shorthand notation for Ci{x,x') and Iq is a Bessel function. 



3 Statistics of the Gabor coefficients 

3.1 Second order statistics 

Second order statistics of the Gabor coefficients of a speckle pattern were calculated in [T7]. 
In this section we study the higher order statistics. In particular we show that all the odd 
moments are zero and that the fourth moment is dominated by the Gaussian contribution; 
i.e. the probability distribution of a Gabor coefficient is 'almost' Gaussian. 

We will often use shorthand notation G for G{w,k,x). From the fact that Tim ([2]) is an 
odd function in (x — xq) it is easily seen that {G)^ = 0. In [T7] it was shown that the variance 
of G is given by 

a^{w, k, xo) ^ - 27)e-(i"^)'"'^' sinh^w^k'^, (16) 

where 7 = i[l + M'^^'^ / {2w^)]-^ and S 1.29 is a numerical constant. The constant S 
originates from an approximation of the correlation (jlOp by a Gaussian curve. The correlation 
between two different Gabor coefficients was also computed. The correlation is defined as 

(^G{w,k,x)G{w',k',x')'j 

Gc{w,w' ,k,k' ,x,x) := — zr^^ (1''') 

(Jg{w, k, x)(Jg{w' , k', x') 
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and the following result was obtained for w' = w, 

\2-\ 



Co 



exp 



7 {x* — X) 



X 



^Wk-k' cos[7(x' - x) ■ {k' + k)] - e"^'""^-^' cos[7(f' - x) ■ (k' - k)] 
2 \J sinh ^w"^ k'^\l sinh k'"^ 



(18) 



The result ()18p is accurate for small distances \x' — x\. For larger distances, the Gaussian tail 
underestimates the actual correlation. 



3.2 Computation of higher moments 



We present a procedure that allows for the computation of arbitrary moments of G. Substi- 
tution of d?]) into ([2]) gives an expression for G in terms of the random phases, 



G{w,k,xo) 



—^e 2 y agUj- ex.p —[b — a + 2xo ■ [a — b)\ 

° AZ 

d,b 



exp 



(a-b)' 



sinh 



2lT 

Yz 



uP'k ■ {a — b) 



(19) 



Taking the nth power of (jl9p and averaging over the random phases leads to a correlation 
function of the form (jl3p . Each 6j evaluates to For each permutation, the sums 

^ (aj — bj) and X^, (&? — a?) reduce to zero. In this way we obtain the following expression 



((G)") 



w 



9 n 



(20) 



t=i 



It is immediately clear that odd moments vanish, since (j2Up is real-valued only for even n. 

We can simplify ()20p by noting that the product of sinh's is zero for permutations that 
have a fixed point. In other words, all cycles of a permutation have to be longer than 1. Also, 
we note that the a-summation factors into a product of independent sums in accordance with 
the cycle structure of the permutation. For instance, the permutation (231) (564) leads to 
the factorisation (Eaiaaas • " OCE 



040506 



, where the dots indicate an expression depending 
only on the three denoted variables. Furthermore, the outcome of each such factor depends 
only on the length of the cycle, and not on the identity of the summation variables. 

The 4th moment is obtained as follows. Among the 4! permutations of {1,2,3,4} there 
are 3 containing two cycles of length two and 6 containing one cycle of length four. The 3 
permutations with two cycles give rise to a contribution 3(0"^)^, which precisely corresponds 
to the 4th moment of a Gaussian distribution. The non-Gaussian contribution from the 6 
remaining permutations (i.e. the 4th cumulant) is computed in appendices lAl and [Bl In the 
regime <C M it turns out that the ratio of the non-Gaussian part to the Gaussian is 1/64. 
In the regime w ^ 3M this ratio is of order 0{M'^ /w'^). 
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3.3 Estimated entropy of a set of Gabor coefficients for given detector noise 



We use a Gaussian approximation for the distribution function of the Gabor coefficients in 
order to derive an upper bound on the entropy of a set of Gabor coefficients, at a given noise 
level of the detector. This is a useful exercise for two reasons. First, it is not a priori clear 
how much of the information present in the source ends up in the Gabor coefficients. Second, 
detector noise affects the Gabor coefficients in a nontrivial way. 

As mentioned in Section [TT2l a relatively small set of coefficients can capture almost all the 
information available in a speckle image. We will consider the example given in Section 11.21 
which is also the choice made in |ll tll7l [8]. Taking more than one width w, or more than one 
wave number \k\ does not make much sense, since all the features in a speckle pattern have 
more or less the same length scale, namely the average speckle size. It is also not very useful 
to take more than two angles of the wave vector: Eq. (jlSp shows us that there is a strong 
correlation sinh[7i(;^A;^ cos C]/ sinh7tt;^A;^ between Gabor coefficients at the same position x, 
when their fe-vectors have a mutual angle C. Hence, for small C there is a lot of redundancy. 

We take a single width w, a single wave vector length k, two perpendicular directions ^i, 
ip2i and a lattice of positions x. We introduce the following definitions: 



Y^Q \x, x') 



G{w, kj,x)G{w, kj,x') 



^avd -7)exp 



7 {x' — x)^ 



2 w"^ 
G{w, ki,x)G{'w, k2,x') 

-/f,(l-27)e(^-i)-''=%xp 2 
X sin[7A;i • {x' — x)] sin[7A;2 • (x' — x)]. 



cos27fc, • (x' - x) - e-^"'''] 



7 (f — x)^ 



(21) 



We define the combined covariance matrix Tiq as 



^[12] 



^G 

V[22] 

^G 



(22) 



We consider the joint probability distribution for all the Gabor coefficients to be Gaussian. 
This, of course, is not true, as we see from the nonzero even moments in Section r3.2[ However, 
for a given mean and covariance matrix, the Gaussian distribution has a higher entropy than 
any other distribution. Hence our procedure yields an upper bound on the entropy. 

We furthermore assume that the detector noise is Gaussian, independent of the intensity 
and independent for each pixel. We can then apply the well known channel capacity formula 
(see e.g. [2]), which expresses the mutual information between a source and a noisy detection 
as the logarithm of a signal to noise ratio. Let G be a vector of Gabor coefficients, and 5G 
the detector noise in these coefficients, then 



I(G; G + 5G) = \ log2 | det(l + S^^Sg; 



(23) 



Here the matrix S^v is the covariance matrix of the noise. We determine S^v as follows. The 
independence of the noise in each pixel yields 



{51{x) 61{x'))n = N]t (5(x - f 



(24) 
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where Nj denotes the noise amphtude, t is the cireci of Si detector pixel and the notation (*)n 
denotes a noise average. Note that Nj is lower bounded by the shot noise. 

Using (I24p and the definition of the Gabor coefficients ([I]) we obtain the following expres- 
sion for the covariance of the noise in the Gabor coefficients 



6G{w,ki,xi) 6G{w,k2,x2] 



Nft / d'^xTiMiw,ki,xi,x)TiMiw,k2,X2,x). 



(25) 



Evaluation of the integral is straightforward (it is equivalent to the derivation of in [T7] 
in the limit S — > 0) and yields, in the special case |A;i| = \k2\ = k, 



^Tv'^ (x, X ) 



^[12] 



{x,x) 



Nft 



Nft (g'-g) 
re 



cos kj ■ {x' — x) 



(26) 



imp- e 2 



sin[iA;i • {x — x)] sin[iA;2 • {x — x)]. (27) 



Sat has the same block structure as Sc. Notice that both Sg and Stv depend on x and x' 
only through the difference {x' — x) . This allows for efficient computation of the determinant 
in (I23p by diagonalisation in the Fourier domain. Notice further that the off-diagonal blocks 
{ki _L k2) have very small values compared to the diagonal blocks due to the presence of the 
sine factors. 

We compute the Fourier transforms as follows. Strictly speaking, the transform is a 
summation over the finite x-grid. We denote the size of the grid as L. However, we will sum 
to infinity, since the error introduced in this way is only an edge effect. The finiteness of L 
is still reflected by the discretisation of the momentum p conjugate to x' — x. All momenta 
are multiples of vr/L. The highest momentum is determined by the lattice constant £ of 
the X- grid. As a second approximation, we will replace the summations over x and x' by 
integrations, i.e. — > / d^x. This is a good approximation provided that the lattice 
constant is significantly smaller than the average speckle size. Using this procedure we obtain 



w 



2 1 



d2(f 



e '^1 sinh w kj • p 

ip-{x'—x) 



(28) 



(29) 



Substitution of (j28p and (j29p into (j23p gives us the mutual information for one Gabor direction. 



+5^^]) = llog2|det('l + (S 



lX]^°S2 (l + cie 



-C2P' 



kr log2 



1 + (S 



N 



^G 



(30) 



where we have defined the constants ci, C2 as 



ci 



' Nft 



27rS 



2 ^av 



Nf t 



C2 



(31) 
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Note that C2 is proportional to the average speckle area, while ci plays the role of a signal to 
noise ratio (SNR). Note too that all dependence on w and k has disappeared from ([30]) . The 
reason is that we are computing a generic upper bound. Computation of the actual amount 
of extracted information will in general depend on w and k. 

The summation domain of is given by {pxiVy) = ihj)x^ with i,j € Z and \px\ < ir/i, 
\py\ < tt/£. Eq. (j30p can be further evaluated by approximating the momentum sum by an 
integration. Notice that the summand in (j30l) only depends on the length of p. Hence the 
computation is simplified in polar coordinates, 



E/(^') dp./ dpyf{p') = - dp 



(32) 



Applying the approximation (f32|) to ([30]) we obtain 



L2 

27r 



C2 



-Dilog(- 



exp C2P 

Cl 



) - |c2p'^ Inci 



P=tt/ L 



(33) 



Here Dilog is the dilogarithm function. Note that (|33p is not expressed in bits but in natural 
units ('nats'), i.e. using the natural logarithm In instead of the base-2 log2. 

Eq. (j33|) can be further evaluated if ci is very large (we call this the 'large SNR' case), or 
when Cl is very small ('small SNR'). We define a quantity y as 



y = cie ' . 
The crossover between the two regimes lies around y = 1. 



(34) 



Large SNR 

In this case we have y ^ 1. Applying the asymptotic relation 

Dilog(x) = x + Oix"^) for |x| « 1 

to (l33l) we obtain 



mci ' 



2£2 



+ 0(y- 



£2^2 ,2 



Cl 



L2' 



(35) 



(36) 



Small SNR 

In this case we have y <C 1. Using (j35p and the asymptotic behaviour 

Dilog(-x) = -7rV6-^(lnx)2 + C'(x~^) for x ^ cx) 

we get 



1 



.M' 



2 n 



(lnci)2 + - + 0{y) + 0{-) + Oijj 



Cl 



(37) 



(38) 



Eq. (j33p is plotted in Fig. [T] on a logarithmic scale; The factor was replaced by vrM^ 
to obtain a result per average speckle area. We see the transition from parabolic behaviour 
to linear behaviour around lav/A'/ = 1. The linear part of the curve is the usual 'log SNR' 
dependence, but here it occurs with a nonzero offset. 
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{q-1 iqO -lqI -lq2 Ni 

Figure 1: Mutual information, per average speckle area, between a noiseless and a noisy vector 
G^^^ of Gabor coefficients, plotted as a function of lav/Nj. The other parameters are: L = 800 
pixels, M = 3 pixels, i = 5 pixels, t = 1 pixel"^. The radius of an average speckle has been 
set, somewhat arbitrarily, to M . 

4 Perturbing speckle patterns 
4.1 Phase perturbations 

We introduce a method for perturbing a speckle pattern in the model of Section [2.11 We shift 
all the phases by random amounts, 

^ ^d = + £d- (39) 

This is done independently for all locations a. The perturbations e are chosen from a dis- 
tribution that is uniform on the interval {—q, q] and zero elsewhere. The parameter q is the 
'strength' of the perturbation. The maximum value is q = n, resulting in a completely inde- 
pendent speckle pattern. The minimum value is q = Aip/2 (see Section [2?2]) . in accordance 
with the uncertainty relation, and represents no visible change at all. The exact relation 
between q and the actual physical perturbation is hard to define. We will come back to this 
in Section O 

We use 'hat' notation (") for perturbed quantities, e.g. I{x) is the intensity after pertur- 
bation. Note that the statistical properties of are exactly the same as the properties of 
ag, i.e. uniformly distributed over the unit circle in the complex plane. Hence the statistics 
of / and G are the same as for the unperturbed quantities. This is precisely what we want 
from our model, for there must be no 'preference' in the formalism for either the perturbed 
or unperturbed speckle pattern. 

Definition ([39]) allows us to study perturbations quantitatively, on a continuous scale as 
a function of q. Of special interest are the correlations between states before and after a 
perturbation. 

We introduce the notation (•)^ for averaging with respect to the random variables {e^}. 
We list several properties that will be useful in later sections. The average effect on is 
multiplication by the following factor: 

(expi.,),= rd.2!^ = 2^. (40) 
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From this the correlation between perturbed and unperturbed a fohows, 

C{as, ag) := ("sag)^^^ = (41) 

Here (•)^ ^ indicates averaging first over e and then over ip. Other useful identities are 

sinz^g 

(expzi/e^) = (42) 

vq 

and 

(exp ^{e, - e^)) ^ = 5, g + (1 - S.J ' . (43) 

For convenience later on we introduce the following shorthand notation, 

0-=^. (44) 

Note that Q = 1 in the case of zero perturbation, and (5 = when the perturbation has 
maximum strength [q = tt). 

4.2 Mutual information in the source plane 

The mutual information [2] between two stochastic variables X and Y is denoted as \{X;Y). 
It represents the amount of overlap in the information they carry, and it can be computed 

I(X; Y) = H{X) + H{Y) - H{X, Y) = H{Y) - H{Y\X). (45) 

This quantity is of great importance in cryptography. Let X and Y be two noisy versions of 
the same secret, possessed by two parties respectively. These parties can derive a common 
secret key from the variable that they hold. The theoretical maximum length of their common 
key is precisely given by the mutual information \{X\ Y). 

In our case, X is the unperturbed speckle source a and Y is the perturbed source a. 
Their mutual information gives an absolute physical upper bound on the length of the key 
that can be derived from a speckle pattern in a reproducible way, for a given noise level q. 
The entropy H[a] is given by ([6]). The conditional entropy -fr[a|a] equals the uncertainty in 
the set {ag} if {as} is known. This is precisely the amount of information contained in the 
perturbation {eg}. Thus we have 



2q_ 

Aip' 



H[a\a] = H[e] = N,,^ log (46) 



Substitution of ^ and (conditional) into (|45]) yields the result, 

I(a;a) = A^reglog^. (47) 

Clearly, if q has its maximum strength (vr) then the mutual information is zero. When q has 
its minimum value, A(p/2, then the mutual information is equal to the full entropy ([6]) of the 
source. 



^ H{X,Y) stands for the combined entropy of X and Y. The notation H{Y\X) denotes the uncertainty in 
Y given knowledge of X. 
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Figure 2: Mutual information between unperturbed and perturbed source, as a function of the 
perturbation strength q. 



4.3 Effect of perturbations on the intensity 

The effect of a perturbation on the intensity is computed as follows. First we average the 
perturbed intensity over the perturbations eg. Using the representation d?]) and the 
identity (j43]) we get 

'i{x)) = + (1 - Q)/av. (48) 



Eq. (|48p shows how, on average, the perturbed value /(x) gradually changes from I{x) to the 
average intensity lav as a function of q. 

Next we compute the two-point correlation function between the perturbed and unper- 
turbed speckle pattern. Taking the expectation value (•)^ of (j48]) multiplied with /(x') we 
obtain 

(i{x)i{x')) -il 

C(/(x),/(x')) := ^ = QCi{x,x'), (49) 

air 



where we have used the fact that the correlation with the constant number lav in the last 
term of (j48|) vanishes. The result (j49|) is intuitive as it states that the total correlation is the 
product of the ordinary two-point correlation Cj ([10]) and a perturbation effect. 
The joint probability distribution p(I, I') is obtained by computing the moments 
/'"(/')"^ \ £qj, ^ gj^^ This exercise is completely analogous to the derivation of (jl4p . 

but now with extra phase factors expieg. Without showing the derivation, we mention that 
the averageing procedures (•)^ and (•)^ in this computation commute and give the result of 
the (^-average. 



■i{x)ni{x')Y 



/"+"^n!m! 2^1 



-n, 



-m; 1; 



iVr, 



. a-{x — x') 



(50) 



Taking the e-average of (|50p is higly nontrivial in general. However, an estimate with errors 
of order l/N^^g <C 1 is easily obtained by replacing averages of powers by powers of averages. 



/!l+"^n!m! 2^1 



-n, 



1 



m;l;QCKx,f) •{l + 0(-— )}. (51) 



reg 
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This gives, up to errors of order 0{N^g^), a joint probability distribution p{I{x), I{x')) of the 
form ([15]) with C/ replaced by QCj. 

4.4 Effect of perturbations on the Gabor coefficients 

The properties of the perturbed Gabor coefficients are readily computed. We define the 
correlation between a perturbed coefficient G{s' , k' , x') and the unperturbed G{s,k,x) as 

G{G, G') := 2 ^— ^, (52) 

with ac as defined in (jl6p . We have again used shorthand notation G for G{s,k,x). The 
averages in ([52]) are easily evaluated. First we note that {G)^ = (^G^ = 0. Second, we 

factorize ^GG"^ = ^G(|g'^ ^ . It directly follws from the definition of G ([T]) and the 
perturbation- averaged intensity (|48p that 

(g')^ = QG'. (53) 

Thus we obtain 

G(G, G') = Q- Gg{w, w', k, k' , x, x'), (54) 

with Gg as defined in (|17p . Hence, as in the case of intensity correlations, the correlation 
between the unperturbed G at location x and the perturbed G' at x' factorizes into a contri- 
bution from the perturbation and the ordinary correlation function Gq- 



4.4.1 Mutual information between perturbed and unperturbed Gabor coeffi- 
cients 

We estimate the mutual information between an n-dimensional vector G of unperturbed 
Gabor coefficients and the vector G of corresponding perturbed coefficients. The vector 
consists of the same set as in Section 13.3] i.e. a single width w, two perpendicular wave 
vectors ki, k2 of equal length, and a grid of points x. We approximate the joint probability 
distribution of G and G by a Gaussian distribution. This is motivated by the results of 
Section 13.21 

We make use of a well known result from information theory. If X is an n-component 
Gaussian-distributed vector with covariance matrix S^, then the differential entropy [2] of X 
is given by 

/i(X) = ilog2[(27re)"|detSx|]. (55) 

Let y be a second vector of the same length, with covariance matrix Sy , and let T,xy be the 
covariance of X and Y. Let Tix and 'Sxy commute. Then the mutual information between 
X and Y follows from the definition (|45p and (|55p. 

I{X;Y) = -llog2 |det(l - S^^S^iS^.^)! . (56) 

Note that in the special case Y = X + N , where N is Gaussian noise uncorrelated to X, ([56p 
reduces to the form ([23]) . 
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Now we consider the case where X = G and Y = G + 6G, where 6G again is the Gaussian 
detector noise discussed in Section [3.31 Note that G has the same (/9-ensemble statistics as G. 
Thus we have 

Sx ^ Sg, Sy^Hc + SAT, T,xY ^ Q'^G- (57) 

Substitution into (j56|) yields the mutual information 

G + 6G) = ^ log2 |det (l + S^^Sg) | - i log^ |det (l + [1 - Q^]^],^^g) \ • (58) 

This represents an absolute upper bound on the information that can be reproducibly extracted 
from a speckle pattern, given that there is noise 61 in the detector and perturbation noise 
{e^} in the source. Clearly, ([55]) reduces to ([231) the case q i 0, and the mutual information 
goes to zero in the case q ^ n. 

We estimate (|58p using the approximation method given in (j32p . This yields, for one 
direction of the fe- vector (and expressed in natural units instead of bits), 



2-KC2 



(1 



1 



-Dilog- 



pC2P 



Cl 



+ Dilog 



(1 - Q2)ci 



(59) 



- p=-k/L 



where the constants ci and C2 are defined in ([ST]) . The result ([Ml) is plotted in Fig. [31 As in 
Fig. [H the factor was replaced by vrM^ to obtain a result per average speckle area. The 
mutual information decreases sharply as a function of the perturbation strength q. 

H (bits 
20 




7T 



Figure 3: Mutual information, per average speckle area, between a noiseless vector G'-'^ of 
Gabor coefficients and a noisy & perturbed version 

Gbl + ^Gbl, plotted as a function of the 
perturbation strength q. The parameters are: L = 800 pixels, M = 3 pixels, i = 5 pixels, 
t = 1 pixel"^ , I^v/^i = 10. The radius of an average speckle has been set, somewhat arbitrarily, 
to M. 



4.4.2 Bit error probability 

We estimate the probability of a bit error in a binarized Gabor coefficient due to the ran- 
dom perturbation. Here we will neglect the detector noise. We use the shorthand notation 
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G = G{w,k,x) and G = G{w,k,x). As in the previous section, we approximate the joint 
probabiUty distribution p of G and G by a Gaussian, making use of the correlation ([5l 



A = a% 



(60) 



1 Q 
Q 1 



A bit error occurs when the signs of G and G are not equal, while |G| > T (see Section [1.2p . 
The probability of this event is given by the following integral expression 

P... = Prob[G < 0|G > n = p,„„G ^ T] = ^G J^dG 

Evaluation of several of the integrals gives 

r-dG ^e-5^'/^G . iErfc^=^ 



-Erfc ^ ' 



(62) 



where Erfc stands for the complementary error function. Fig. |4] shows the behaviour of Perr as 
a function of q and T. Exact evaluation of the leftover integral in (j62p is difficult in general. 
However, in some limiting cases analytic results can be obtained. For instance, for T = 
the result is Pprr = vr^^ arccos Q. Furthermore, in the two limiting cases , ^ 

and Qi approximations can be obtained. In the former case we apply a large- 

argument asymptotic expansion of the Erfc function; in the latter case a Taylor expansion. 
For the first case, let us define the small parameter e ^ 1, 



After some straightforward but tedious algebra we obtain 

P^^^ = £ ! {e^ _ £5(3/2 + 1tV4) + Oie')} . (64) 

Erfc[T/((TG\/2)] 27raG^ ^ 2 / g; V ;/ v ; 

Expression (|64|) is useful in the weak perturbation limit Q ^ 1 and in the limit of large 
thresholds, T/ac ^ 1- 

For the second case we write r] := 1/e <^ 1. A Taylor expansion of the Erfc function in 
(f62|) yields, after some algebra. 



p — _ 



3 



2fTGi.crG\/2 e 2 



(65) 



Expression (|65|) is useful in the strong perturbation limit Q ^ and in the limit of small 
thresholds, T/ac ^ 1- 
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Figure 4: Bit error probability Pcrr according to \6S^) as a function of the perturbation 
strength q for various values of the threshold T. 

5 Comparison to experimental results 

In this section we briefly compare a number of theoretical results, obtained in the previous 
chapters, to actual experiments. The experimental data were obtained with a very simple 
setup, consisting of a laser, a sample holder and a detector. The laser has a wavelength of 
A =780nm. A parallel beam shines on the sample at an angle of 45°. The spot is circular, 
with a diameter D =lmm. The sample is a piece of paper. The detector is a CCD camera, 
mounted at a normal angle to the sample. The distance between the sample and the camera 
is z =10cm. The camera has a pixel pitch of 6.25/um and takes 1024x768 pixel images 
with 256 gray scales. The typical speckle diameter at the location of the camera is of order 
Xz/D =78/um, corresponding to 12 pixels in the image. The setup is not particularly well 
protected against background light. 

Intensity distribution 

In order to illustrate the quality of our data, we show in Fig. [5] a histogram of the gray 
levels present in a single typical CCD image. The lowest gray scale present in the image was 
normalised to zero. The deviations from the theoretical curve (|12p at low intensity show that 
there is a noticeable effect of the background light. 

Statistical distribution of the Gabor coefficients 

Here we discuss the experimental verification of the theoretical results of Sections 13. II and 13.21 
First we show that the theoretical prediction (jl6p for ac is accurate. Then we show that the 
distribution function of the Gabor coefficients has a noticeable deviation from the Gaussian 
form, with fatter tails than a Gaussian. 

The empirical (/3-ensemble probability distribution of G{w,k,XQ) should ideally be ob- 
tained as follows: Insert many samples; for each sample, measure G{w,k,X()); finally make 
a histogram of all the Gabor coefficients; this yields the empirical distribution function for 
G{w,k,xo). 
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Figure 5: Intensity histogram of a single speckle pattern. The histogram was made 
with a bin width of 5 gray values. The theoretical curve is the exponential distribution 



We used a less labour-intensive approach. We took a single sample; from the single CCD 
image we derived Gabor coefficients G{w, k, xq) for all xq; we made a histogram of the resulting 
set; we used this as the empirical distribution oi G{w,k,xo). This approach is motivated by 
(a) the fact that the (/^-ensemble probability distribution does not depend on xq and (b) the 
ergodicity property of laser speckle, i.e. the property that the spatial intensity distribution 
asymptotically tends to the (^-ensemble intensity distribution. 

The result is shown in Fig. [U where ac is given as a function of k = \k\ for a number 
of choices for w. The theoretical result (jl6p is also plotted. The correspondence of theory 
vs. experiment is very good, except at large k. There we start to see the difference between 
the spatial continuum appraoch of the theory and the discrete pixellated nature of the CCD 
images. The theory uses spatial integration, while the data processing involves summation 
over pixels. For fast oscillations, the integration in ([T]) averages out to zero more quickly than 
the summation. 

Fig. \7\ shows the shape of the empirical distribution function of the Gabor coefficients. 
The curves were derived from a single image. It can be seen that the distribution has fatter 
tails than a Gaussian, as was derived in Section [3.21 

Effect of perturbations 

We investigated perturbations as follows. Our sample was a piece of paper whose surface 
structure changed over time. We took 40 pictures at half hour intervals. The changing surface 
structure can be regarded as a random perturbation as modelled in Section [4Tl Unfortunately, 
it is not possible to experimentally regulate the perturbation strength q. We therefore used 
the following approach to compare theory and experiment. We looked at all pairs of CCD 
images; there are {^2) — '^^^ pairs. For each pair (A, B) we computed the empirical intensity 
correlation H/ and Gabor coefficient correlation He, 



(l/4v)exp(-//4v). 



(66) 
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Figure 6: The variance ac of the Gabor coefficients (normalized w.r.t. lavj as a function 
of the spatial frequency parameter k, for several values of the Gaussian width w. (Here w is 
measured in pixels, and k in radians/pixel.) The solid curves represent the experimental data. 
The dotted curves are the theoretical result [T6\). with the parameter choice M = 5 pixels. 




Figure 7: Distribution of Gabor coefficients obtained from a single speckle pattern, by applying 
two perpendicular angle parameters. Left: (p = 0; Right: ip = 90°. The solid curve depicts a 
histogram of the Gabor coefficients. The dotted curve is the least-squares Gaussian fit to that 
histogram. 



ig[A,B\ = . ' I (67) 

^^pix EJGf ]2 - [iv .1 E. Gt]\ N-^Y.AGf?-[K^i:,G^f 



Here A'pix is the number of pixels in the image, Ei stands for summation over all pixels, if 
denotes the intensity in image A at location Xj, and Of is shorthand notation for G{w^ k,Xi) 
for some fixed value of w and k. The empirical correlation should be equivalent to the 
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theoretical correlation (j49p with x' = x; Similarly, Eg should be equivalent to (j54p with the 
substitution w' = w, k' = k and x' = x. Hence, we expect Ef[A,B] = Q and ^(^[A, 5] = Q. 
In Fig. [5] we have plotted vs. H/ for all image pairs. The data points are clearly bunched 
together on a narrow band slightly above the theoretically expected line = H/. We 
hypothesize that this small difference is due to detector noise, which we did not take into 
account here. The Gabor coefficients, resulting from a spatial sum, are less sensitive to 
detector noise than the intensity itself. Hence the correlation of the Gabor coefficients is 
larger than the intensity correlation H/. 
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Figure 8: Effect of a perturbation on the intensities and the Gabor coefficients. Horizontal 
axis: the intensity correlation function 166\) . Vertical axis: The correlation function l^67\ ) of 
the Gabor coefficients. Data from two perpendicular k vectors are plotted together. The solid 
line is the theoretical predicition for zero detector noise. 



6 Summary 

Laser speckle has been proposed in the security literature as a source of high-entropy bit 
strings for various (cryptographic) purposes. It is important to know how the physical prop- 
erties of speckle affect the entropy of the extracted bit strings. More in particular, we need 
to know the mutual information between two repeated key extractions when noise is taken 
into account. Another important parameter is the bit error rate. 

In this paper, we have developed a simple approach to address these issues. We have stud- 
ied the case of key extraction using Gabor coefficients. We used a simple model for speckle, 
generated by a large number of independent random phases in a source plane. We have mod- 
eled perturbations of the speckle pattern as small, uniformly distributed perturbations of the 
random phases. Detector noise was modeled as being Gaussian, independent of the intensity 
and without correlations between the detector's pixels. 

Our main results are 

• The Gabor coefficients have a distribution function that is close to Gaussian. 
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• We have derived an expression for the mutual information between an unperturbed and 
perturbed speckle source. 



• We have obtained analytical expressions that give an upper bound on the mutual entropy 
of a set of Gabor coefficients (i) when there is detector noise but no perturbation and 
(ii) when there is detector noise as well as a perturbation. 

• We have computed the bit error rate caused by perturbations of a speckle pattern. 

Experimental data on the statistics of Gabor coefficients and on the correlation functions of 
Gabor coefficients and intensities are in accordance with theory. 

The results of this paper, particularly the mutual information and error rate expressions, 
provide useful parameters for key extraction systems. 
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A Fourth moment of the Gabor coefficients; case w > M 



In this appendix we calculate the non-Gaussian part of (G^) . As discussed in section [3l 
(113p has 6 permutations with one cycle. These are equivalent because of the possibility of 
relabeling the dummy variables oi, • • • , 04. Introducing the notation Sij 
we can write the non-Gaussian part of (G^)^ as 



smh^k-{ai-aj). 



iX^ 



-2w'^k 



2 1,2 



W 



exp 



w 



[ai • a2 + 02 • 03 -I- as • 04 -I- 04 • aij 



5'l2'S'23'S'345'41- 



(68) 



Note that the expressions in the exponents are invariant under relabeling of the summation 
variables. This allows us to expand the product <S'i2'S'23S'345'4i into sixteen terms, which can 
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then be grouped together into contributions of the same 'type', i.e. equivalent under the 
summation. This gives 



5'l2'S'235'34'S'41 



1 



1 ' 



«2) +cosh^fe • (ai 



^2 + a 3 — 04) 



(69) 



Next we diagonahse the quadratic terms in the exponent of (j68p . The exponent is of the form 
exp(— /3^a/xa''"), where a denotes the four-component row vector (oi, • • • ,04), /? = w/{MR), 
and the matrix /i is given by 
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(70) 



/X has one eigenvalue with eigenvector (1, 1, 1, 1). The other three eigenvectors are 

^(1,0,-1,0) -i=(0,l,0,-l) 1(1,-1,1,-1) (71) 

with eigenvalues 1, 1 and 2 respectively. We defining the new summation variables vq = 
P{ai + a2 + as + 04), vi = /3(ai - as), V2 = /3(a2 - 04) and us = /3(ai - 02 + as - 04). Next 
we approximate the summations by integrals: — > J d^a. Taking the Jacobian into 
account, we write / d^ai • • • d^a4 = / d^uo • • • d^vs. Thus (p8]) can be approximated as 





g-2«>2fc2 / 







d^VQ 



[1 + cosh 2wk ■ V3 — 2 cosh wk 



•••d\s exp [-^{vf + + vD] 

V3 + V1- V2)]- 



(72) 



Given the finite summation intervals of ai---a4, evaluation of the integrals in (j72|) does 
not yield esthetic results. However, if u; ^ 3M then a substantial part of the Gaussian 
distribution is covered by the integration, and considering the interval to be infinite is not a 
bad approximation. The case w < M is discussed separately in Appendix[Bj We only have to 
take into account the finiteness of the vg-integral. The integrand in (|72p does not depend on 
vq, and this leads to a factor 7r(4/3i?)^, since |iTo| < 4/?i?. The remaining integrals are readily 
evaluated. The final result is 



34 + §C 



w 



2w)2fc2 



(73) 



Prom ()16p we see that is asymptotically proportional to l'^-^{M/w)'^ for small M/w. Hence 
the result ([73]) is of the form 3o-^(l + 0[M'^/w'^]). 



B Fourth moment of the Gabor coefficients; case w <^ M 

In this appendix we compute (G^)^ in the limit where the length scale w of the Gabor 
transform is very small compared to the average speckle size. In this limit, the intensity 
changes only slowly as a function of x within the Gaussian envelope. Around the point of 
interest xq we can make a linear approximation 

I{x)^I{xo) + D-{x-xo), (74) 
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with D = V/(xo)- The Gabor transform ([T]) then reduces to 

Differentiating ([7]) and taking the inner product with we obtain 

D ■ k = ^' aga*-k ■ (a — b) exp — iP — + 2xq ■ (a — b) 

y6 b \^ 



d,b 



(75) 



(76) 



We square (j75p . apply Q and replace the sums by integrals. In this way we obtain 



= {G''{w,k,xo)) ^ 4w^k''llM-\-^ 



(77) 



The fourth moment of (j75p is obtained using (jl3p and again replacing summations by inte- 
grations, 



G^(«;,fc,fo)) ^3(7^(1 + ^). 
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